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Abstract -The modeling framework of several recent studies of granular impact is a force law 
consisting of three terms: gravity, a depth-dependent static term, and a collisional term. We recast 
this usual force differential equation for F(t) into an equation for the kinetic energy vs. depth, 
K(z). We show that this reformulation is beneficial, both theoretically and experimentally. With 
the original force-law form, obtaining closed-form solutions may be difficult, and experimental 
comparisons require acceleration data, which is difficult to obtain at high speeds. By contrast, 
with the kinetic energy formulation, standard differential equation methods yield solutions for 
K(z). We apply the kinetic energy approach to experimental data for the trajectory and velocity 
of an intruder that is impacting on a quasi-two-dimensional array of photoelastic disks. We also 
directly relate the two phenomenological terms in the force law to detailed properties of the 
granular medium and intruder. 
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Introduction. — A dense granular material which is 
struck by a high-speed object exerts a decelerating force on 
the intruder, the nature of which is important for many 
applications, such as soil penetration tests [TJH], meteor 
impacts [3], and ballistic applications [HE]. Additionally, 
understanding the dynamics of intruder motion in a gran- 
ular material, as well as the granular flow around it, is a 
fundamental problem in granular physics. To probe this 
process, it is common to drop or push an intruder into 
a granular material and record the dynamics and/or the 
force, F, exerted on the intruder. In general, the dynam- 
ics depend on the microscale physical characteristics of the 
grains and intruder, and may show large fluctuations, as 
is common in granular materials. A complete description 
represents a complex problem in granular rheology. 

A common approach [61414] dating back to the times of 
Eulcr and Poncelet is to use space-time averaged macro- 
scopic force laws, where the various terms in the law are 
empirical expressions based on assumptions of relevant 
physical principles. These terms are typically assumed to 
depend on the depth, z, as well as on the intruder velocity, 
v, taken to be strictly vertical, i.e. v = z. A fit of a given 
model to experimental data provides a way to explore the 
interplay of the granular material and the intruder, and 
the roles of various parameters. Although these models 
are empirical, they are often reasonably successful in de- 
scribing the average dynamics of intruder trajectories. For 



the case when the initial velocity of the intruder is moder- 
ate, a number of these models, have the following general 
form: 

F = mz = mg — f(z) — h(z)z 2 , (1) 

where dots denote time derivatives. This force law con- 
tains three terms: gravity, a depth-dependent static term, 
often taken to be linear in z, and a collisional term pro- 
portional to the square of the intruder speed. Here, z 
is measured from the original unperturbed surface, with 
z = as the point of initial contact. These force laws are 
coarse-grained descriptions of local granular processes in 
a certain physical regime (in this case, moderate but sub- 
sonic intruder velocity), so the specific forms of f(z) and 
h{z) vary, depending on microstructure and preparation. 

Fitting to these models is often done by examining data 
from trajectories with many different initial velocities [lOf - 
IT4"] . For a specific depth, z = £, eq. ([T]) becomes: 



F = mz = mg — /(C) — h(()z 2 



(2) 



If this model is valid, plotting z versus i 2 , calculated at 
z = C, yields an approximately straight line of slope h(() 
and intercept g — f(C)/m. 

Applying these models, experimentally or theoretically, 
can be somewhat complicated. Obtaining closed-form so- 
lutions for the trajectory requires specific assumptions 
about f(z) and h{z). Also, it is often difficult to obtain 
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accurate acceleration data at short time scales, since ac- 
celerometers have limited time-resolution, and discrete dif- 
ferentiation from position data greatly amplifies measure- 
ment noise. Additionally, large, high-frequency, physical 
fluctuations in instantaneous acceleration for the exper- 
iments considered here |14j are absent in the algorithm 
described above. 

In this letter, we suggest a new approach to this model, 
which reformulates eq. ([2]) in a kinetic energy picture, and 
we show that this approach is beneficial in regards to both 
theory and experiment. Mathematically, it allows formal 
closed-form trajectory solutions without specific assump- 
tions on f(z) and h(z), and it requires only data for veloc- 
ity vs. depth, which is easier to obtain. Furthermore, data 
for v(z) has smaller physical fluctuations than acceleration 
data. 

We apply this approach to two-dimensional impact ex- 
periments using photoclastic disks, recorded by a high- 
speed camera (40,000 frames per second), which yields the 
dynamics of the intruder as well as the local granular force 
response. Working within the kinetic energy framework, 
we use the intruder dynamics to determine f(z) and h{z) 
in eq. ([TJ, for our experiments, and hence to infer infor- 
mation on grain-scale mechanisms. 

Kinetic Energy Formulation. — We recast eq. ([T]) 
from second order in time for z to first order in depth 
for kinetic energy of the intruder, K = ^mz 2 , by using a 
relation that is familiar from the work-energy theorem of 
mechanics: mi = dK/dz. 



dK 

dz 



mg - f(z) 



2h{z) 



K. 



(3) 



This yields an inhomogeneous linear ordinary differential 
equation with (potentially) nonconstant coefficients, by 
contrast to the nonlinear equation of motion, eq. [TJ for 
z(t). The fact that eq. [3] is a linear ODE means that stan- 
dard ODE methods immediately yield formal solutions for 
K(z): 

K(z) = K p (z)(K + d>(z)). (4) 
where Kq is the kinetic energy at impact, 



and 



K p (z) = exp I - / (2/m)h(z')dz' 



<P= / dz'[mg-f(z')]/K p (z'). 
Jo 



(5) 



(6) 



This reduces the problem for the trajectory to a quadra- 
ture. The velocity can be written as: 



dz 



-Kp(z){K + cj>(z)) 
m 



z(t), follows by integrating and inverting: 



t(z) = / dz 1 
Jo 



m 



K p (z')(K + <j>(z')) 



1/2 



-1/2 



(7) 



(8) 



If the forms of f(z) and h(z) are simple, much of the 
calculation of these integrals can by done explicitly. For 
example, using the commonly assumed forms f{z) — fo + 
kz and h{z) = b, wc obtain K[z) as 



K(z) = (K - ci) exp(-c 2 z) + Ci 



c 3 z. 

-k]/c 2 \ 



(9) 



C2 



Here, the constants are c\ = [(mg — /o)c 2 
2b/m, and C3 = kjc-2 = km/(2b). 

Even without integrating eq. (JHJ) , it is possible to find the 
stopping distance by setting K(z stop ) — 0, or <p(z s top) = 
—Kq, yielding the stopping depth as a function of impact 
energy, Kq. Specifically, for the common case described 
by eq. ©, the stopping depth, 



'stop 



satisfies 



Zstop 



2b 



log 



26 TS 



{fo + kz s top 



km 
26 



mg 



(10) 



Note that if we take f(z) as roughly constant, f(z) = fo 
and k = 0, then z s t op increases logarithmically with Kq, 
as in [5]. 



*stop 



2b 



log 



1 



2b 



m V fo - 



(11) 



This approximation is relevant in the limit of high im- 
pact energy, where bz 2 dominates. In the low energy 
limit, it predicts z s t op (Ko = 0) = 0. However, an in- 
truder must be at least partially submerged to be sup- 
ported by frictional grains. This occurs at a depth which 
increases with intruder size. So, when Kq — > 0, we expect 
z s top{Ko = 0) ~ D. Note that eqs. ^ and ([TUf may yield 
a nonzero upward force, F = ^p, as the intruder comes 
to a stop, and we return later to this issue. 

Measuring f{z) and h(z). This formulation also pro- 
vides a way to measure f(z) and h(z) in terms of ex- 
perimental data for z(t) and i(t), without specific as- 
sumptions for f(z) and h(z). Subtracting two different 
trajectories, with different Kq's (not necessarily close), 
Ki(z) = K p (K i0 + 4>) and K,(z) = K p (K j0 + 0), yields: 



Ki(z) 



K i0 
which gives, 



z •2h{z')dz' 



K 



jo 



h(z) = - 



d \m 



log K p {z) 



(12) 



(13) 



To examine average behavior and scaling, we typically use 
J h(z)dz in our discussion below, to avoid numerical dif- 
ferentiation. 

The expression for f(z) follows from setting eq. [4] equal 
to zero at z = z s t Q p, then differentiating with respect to 
Zstop, which yields: 



stop ) 



dK Q 

dz. 



stop 



mg. 



(14) 



This analysis requires a determination of K p (i.e. h(z) is 
determined), and ;z s t op (-Ko), where the latter is generally 
straightforward. 
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Application to Experimental Data. To test 
the approach outlined above, we use data from two- 
dimensional granular impact experiments |14j , where disk 
and ogive intruders (i.e. intruders with an elliptical-like 
nose), cut from bronze sheet, are normally incident on a 
collection of approximately 25,000 bidispcrse, hard, pho- 
toelastic disks (diameters of 6 mm and 4.3 mm) confined 
between two 0.91 m x 1.22 m x 1.25 cm acrylic sheets. 
Particles are constructed from PSM-1, a stiff photoclastic 
polymer, made by Vishay Measurements Group. Intruders 
have initial velocities between and 6 m/s, which are in 
the subsonic regime, since videos show the granular sound 
speed to be about 300 m/s. We record results with a 
high-speed camera, typically at 40,000 frames per second, 
which allows us to determine the position of the intruder 
in each frame, as well as the forces on the photoelastic 
particles. Here, we focus on the intruder dynamics only; 
sample trajectories are shown in fig. [TJ 

Once the intruder trajectory is known, we differentiate 
to find the velocity and acceleration of the intruder at 
each frame. To avoid noise amplification, some filtering 
is required with each derivative. This is accomplished by 
fitting a linear function to W frames of the position data, 
centered at z(t) for the frame of interest, which yields the 
velocity, with time resolution reduced by a factor of W. 
The same procedure is repeated to obtain the acceleration 
from velocity data. We choose W = 300, which is the 
smallest value for W such that the signal-to-noise ratio is 
10:1 for acceleration data. This ensures that the observed 
fluctuations in velocity and acceleration are physical. As 
discussed in Clark et al. [14] , the acceleration data for the 
intruder, obtained in this manner, exhibits fluctuations 
that are intrinsic to the emission of acoustic pulses at the 
interface between the intruder and grains. 

Here we use eight different intruders, with width D and 
nose shape C(x), shown in fig. [TJ Four circular intruders, 
with diameters, D, of 6.35, 10.16, 12.70, and 20.32 cm, 
were used to test size effects, and four ogive intruders were 
constructed to test shape effects. The shapes of the ogives 
consisted of a continuous piece of material, where the lead- 
ing part is a half-ellipse truncated along the minor axis, 
with semi-major axis a and semi-minor axis 6 = D/2, ter- 
minated by a rectangular tail of length L. Three different 
ellipses were used, with a/b = 1 (half-circle), a/b = 2, 
and a/b = 3. The width of the ogives was held constant, 
D = 9.3 cm, and L was varied to keep the intruder mass 
constant (L = b = 4.15 cm for a/b — 3 case, longer for 
other ogives). By keeping the width and mass constant, 
we isolate shape effects. Additionally, we used one smaller 
ogive with a/b =1,6 = 3 cm, and L = 7.7 cm. 

As discussed above, fitting to the force law of cq. flT]) 
requires plotting the acceleration versus velocity squared. 
As is obvious from fig. [TJ the fluctuations in acceleration 
are large. Thus, determining a clear value for f(z) and 
h(z) at each depth is difficult. However, with the kinetic 
energy approach, only velocity data is needed to deter- 
mine f(z) and h(z). Using eq. (fT2|). and averaging over all 
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Fig. 1: (TOP) Single trajectories (depth, velocity, and accel- 
eration) for the three smallest circular intruders with similar 
initial impact velocities, where t — is the first contact with the 
granular surface, measured from photoelastic response. (BOT- 
TOM) Shapes of all intruders, drawn to scale, as described in 
the text. The close up of the a/b = 3 intruder shows a sketch 
of a simplified single particle collision, which can be used to de- 
rive a form for the scaling behavior of h(z), as discussed later 
in the text. 

pairs of trajectories (omitting trajectory pairs with very 
similar initial velocities), we obtain a clear average value 
for K p (z), and thus for h(z), as shown in fig. [U 

Data for — ^ log K p = J h(z)dz are approximately lin- 
ear in z. The slope gives the collisional term h(z), which 
scales approximately with the intruder size, as discussed 
below. We observe some variation in h(z) with z: an ini- 
tial transient regime, after which h(z) approaches a nearly 
constant value. We note that for circles, the collisional 
term is stronger right at impact, which may be surprising, 
since the area of impact is smallest then. This effect is 
reminiscent of surface tension or so-called "added mass" 
effects for fluid impacts [T8H2T)] , but it is not obvious what 
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Fig. 2: (TOP) Plot of / h(z)dz = -{m/2)\ogK p for all eight 
intruders. All are approximately linear, and the local slope 
gives the value of h(z). Inset shows K p (z) = AK(z)/AKo, 
computed for all pairs of trajectories, i and j, with a minimum 
AK = K i0 - K jQ . (BOTTOM) Taking the local derivative 
shows transient behavior, which depends on the shape of the 
intruder. Circular intruders show stronger h(z) near the sur- 
face, reminiscent of similar effects in fluid impacts |18H20j . For 
the ogive shown, this effect is reversed: h(z) starts lower and 
increases sharply to a local maximum, then weakly decreases. 



physical mechanisms are at play in the granular case. 

Once h(z) is specified, cq. (fl~4| provides a solution for 
f(z) with adequate data for z stop {Ko). First, we plot the 



final depth, 



'Stop; 



versus An, as shown in fig. |3l Note 



that the data for the higher energies are consistent with 
the logarithmic behavior in eq. (fTTj) . and that the data 
for low energies deviate from this curve with a non-zero 
intercept which scales with intruder diameter (z stop (A'o = 
0) « D). Only circular intruders are used, since ogive 
intruders penetrate much deeper and may interact with 
the lower boundary of the experiment. 

Finally, by fitting a smooth function to the curve for 
each circular intruder and differentiating, we solve for 
f(z), wherever z stop (Ko) is defined. This yields a lin- 
ear function for all circular intruders with a non-zero in- 
tercept, /o, which scales linearly with the intruder mass 

(fig. SJ). 

As discussed previously, the net force must go to zero 
as the intruder stops. With this in mind, we examine 
trajectories for v < 0.3 m/s (shown in fig. |4j, where 
h(z)v 2 <C mg. This shows approximately constant de- 
celeration as the intruder comes to a stop (consistent with 
the expected value of f(z) from fig. [3]), and the accelera- 
tion jumps to zero at v = 0, as in [TTIH^] . 
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Fig. 3: (TOP) Plot of z stop vs. K , with fits of the form 
alog(&A + 1) + c. (BOTTOM) Plot of f(z) for circular in- 
truders. Linear fits are fo + kz, where the slope, k, corresponds 
to hydrostatic pressure. Note that f(z) is dominated by the 
offset, fo, for our data. Inset shows plot of fo vs. mg, with a 
linear fit through the origin, with slope 1.35. 
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Fig. 4: Main plot shows the average velocity of all four cir- 
cular intruders as they come to a stop. All intruders decel- 
erate at the same rate, a w 0.75g. The upward force re- 
quired would be approximately 1.75mg, which is consistent 
with f(z 3 top) = 1.35mg + kz a t op , from fig. [3] Inset shows the 
end of all trajectories, with initial velocities between 1 and 
6 m/s, for the D = 12.7 cm circular intruder. Note that the 
stopping time is very weakly dependent on initial velocity, since 
these trajectories all end at approximately the same time. 
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Connecting to microstructure effects. The 

terms f(z) and h(z) are intended as coarse-grained de- 
scriptions of local granular processes, similar in spirit to 
the static and collisional terms in a general coarse-graining 
description, such as that formulated by Goldcnbcrg and 
Goldhirsch [3T]. With this in mind, we now examine the 
behavior of f(z) and h(z) with the goal of understanding 
the specific kinds of grain-scale mechanisms in our system. 

First, for h(z), photoclastic analysis shows that, for 
moderate to high velocities, the intruder deceleration is 
dominated by collisions with grain clusters, which send 
energy away in the form of acoustic pulses along networks 
of grains [l4j . We expect that the mass associated with 
grains in the collision is greater than just a single grain, 
since a collision involves forcc-chain-like structures that re- 
main in contact over a finite time. A similar collisional pic- 
ture was recently proposed by Takehara et al [T2J. Thus, 
we consider the simplified case shown in the bottom right 
corner of fig. [U where an intruder with mass m, width D, 
and velocity v collides with particles (or particle clusters) 
which have mass m g and are at rest. We assume energy 
loss during each collision which is captured by a restitu- 
tion coefficient, e. In this case, each collision imparts a 

-v cos 9, where n is 



momentum change An = h(l + e)- 

° r \ j rn g -\-m 

the unit normal to the intruder surface and 9 is the an- 
gle between the velocity and n. We take the typical time 
between collisions as At = ad/(v cos 9), where d is the 
particle diameter and a is an 0(1) constant. Thus, the 
average force in the n-direction is: 



/ = 



Ap (1 + e)v 2 cos 2 



At 



ad 



(15) 



The net force is given by the number of these collisions, 
which in a length dl is given by dn = /3dl/d, so dF = 
hf dn = (hf/3/d)dl. If the shape of the intruder is given 
by z = C(x), as in fig. [U then dl = (1 + C*' 2 ) 1 / 2 , and 
cos 6* = (1 + C" 2 )" 1 / 2 . We integrate over the leading edge 
of the intruder, keeping only the z-componcnt of F: 



F z = 



dF ■ z = 



= B TI[C(x)] ■ v 2 - = h{z)v A (16) 

(1 + e)ir/3/(ad 2 ), T = m g m/(m g + m) is the 
reduced mass, and the shape factor I[C(x)} is given by: 

r-D/2 



Here, B 



fcos9(P/d)dl 
- 2 = h(z)v 2 
- m g m/(m g 



I[C{x)] 



D/2 



dx{\ + c u y 



(17) 



Given C'(x), we can evaluate this integral to see how 
changes in shape should affect this force. Note that for 
a given shape (i.e. circular nose), I[C(a;)] scales linearly 
with D, so the width-independent shape factor is given by 
J = I[C(x)]/D. For ogives with circular/elliptical noses, 
the form of C{x) is given by: 



C(x) 



1/2 



(18) 




Fig. 5: Rescaled plot of / h(z)dz/(J(s)D 2 



z/D for all 



intruders. Normalizing by the shape effect, J(s), and the 
intruder width, D, captures the difference between elliptical 
noses and circular noses of different sizes. The two heaviest 
intruders with circular noses, with m/D ~ 3 kg/m, do not col- 
lapse as well. A suggested explanation for this effect is in the 
text. 



where s = a/b is the ellipse aspect ratio. Inserting this 
form into eq. (j!7p and evaluating, we obtain: 



J(s) 



IJs) 
D 



1 



s 2 arctan ys 2 



1 



(19) 



If we rescale data for J h(z)dz by I(s) = J(s)D (fig. O, 
we observe that all intruders except two collapse nicely 
onto a master curve, with BqT ~ 7. The two intruders 
that do not collapse as well are the larger, circular ogive 
and the largest disk. We believe that the deviation of these 
two intruders relates to their larger ratio of m/D ~ 3 kg/m 
(for all other circular- nosed intruders, m/D < 2 kg/m). 
These two intruders generate photoclastic activity that 
extends considerably farther into the material (the larger, 
elliptical- nosed intruders with m/D ~ 3 kg/m generate 
far less photoclastic activity, perhaps keeping them in the 
same regime as the smaller circles). Thus, as a potential 
explanation, we examine the effect of the mass of a typical 
grain cluster involved in a single collision: 



mD 



B J(s) 



1 



1 + m/m £ 



B J- 



(20) 



Here, we assumed m 3> m gi since the bulk mass density 
of intruders is seven times that of particles, and the area 
of particles involved is presumably less than that of the 



intruder area. Thus, we see h/(JD) 



and fig. [5] 



shows that h/(JD) collapses for all intruders except the 
two in question, where it approximately 1.3 times bigger. 
This would require m g to be about 30% bigger for these 
two intruders, which is consistent with observations from 
high-speed photoclastic videos. We note that other ex- 
planations are possible, such as changes in restitutional 
losses or effects from the container boundary. Similar col- 
lective effects were used in [13], as well as suggested by 
Waitukaitis and Jaeger to explain the more extreme case 
of shear thickening of suspensions which are subjected to 
impact [22]. 
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Fig. 6: Plot of rescaled f(z) for all circular intruders, in di- 
mensionless units, with f = 'ymg subtracted off. The fit line 
shows 4>p g , the slope for hydrostatic pressure. 

Next, we consider f{z), which captures depth- 
dependent effects which arc independent of velocity. From 
fig. El we see that f(z) is well captured by a linear function 
with a slope which is comparable to hydrostatic pressure, 
4>p g gz, and an intercept which increases with the intruder 
mass (here, <j> is the packing fraction of the granular mate- 
rial, assumed to be 0.8 in our system). The nature of the 
offset term, f , is unclear, but can be thought of as a bulk 
material stiffness, f = E D, and the linear increase cor- 
responds to stiffening of the granular material at greater 
depths. Thus, we propose a form for f(z) for our results: 

f(z)=E(z)-D = (Eo + <t> Pg gz)D, (21) 

with Eq = jmg/D, where 7 ~ 1 .35. (We note that impact 
experiments performed by Goldman and Umbanhowar |f lj 
show a similar result for f , but with a 7 which increases 
from to 2 during < z < D, then saturates at 7 ~ 2.) 
Fig. [6] shows dimensionlcss f(z) for all circular intruders, 
with offset fo = jmg subtracted. 

Conclusion. — We have shown a new approach to 
a commonly used model for describing the dynamics of 
an intruder impacting on a granular material. By refor- 
mulating the model into a linear ODE, we obtain formal 
solutions of the trajectory in terms of the initial kinetic 
energy, as well as a systematic way of calculating the dif- 
ferent terms in the model (namely f(z) and h(z)) with 
only position and velocity data, which is more easily ob- 
tained experimentally than data for acceleration. 

Additionally, we have used this approach to measure 
f(z) and h(z) for experimental data. Photoelastic analy- 
sis [H] suggests a collisional/acoustic picture, and we have 
shown that a straightforward collisional calculation gives 
a good prediction for how intruder size and shape affect 
h(z) for our experiments. We also observe that the usual 
assumption that h(z) is constant applies only after an ini- 
tial transient at impact, which varies with intruder shape, 
and could be important in engineering and control appli- 
cations. We also note that terms linear in velocity have 
been proposed in the context of this model QT|, but our 
data shows no need to implement them here. Future work 



may include investigation of other intruder shapes (e.g. 
triangular/conical noses), as well as connection to other 
regimes, such as the slow drag regime. 
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